Exotic quantum spin models in spin-orbit-coupled Mott insulators 
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We study cold atoms in an optical lattice with synthetic spin-orbit coupling in the Mott-insulator 
regime. We calculate the parameters of the corresponding tight-binding model using Peierls sub- 
stitution and "localized Wannier states method" and derive the low-energy spin Hamiltonian for 
fermions and bosons. The spin Hamiltonian is a combination of Heisenberg model, quantum com- 
pass model and Dzyaloshinskii-Moriya interactions and it has a rich classical phase diagram with 
collinear, spiral and vortex phases. 

PACS numbers: 



Since the first experimental realization of Bosc- 
Einstcin condensate (BEC), cold atoms have proven 
to be an excellent playground for studying many-body 
physics 0, Q and many interesting phenomena take place 
for strongly interacting atoms in an optical lattice. These 
studies started with an experimental observation of su- 
perfluid to Mott-insulator phase transition Q which was 
followed by experimental and theoretical work on both 
Bose and Fermi gases in lattices of different dimension- 
ality and in various parameter regimes 0,0. The key 
features of cold atoms in optical lattices are the excellent 
tunability of parameters and the fact that the sample is 
almost perfectly described by the Hubbard model in the 
deep lattice regime Q. Since it is well known that the 
Hubbard model is mapped to an effective spin Hamilto- 
nian for Mott-insulator phases with integer filling it 
is clear that cold atoms can be used to "engineer" various 
quantum spin systems [H], Q from those described by the 
Heisenberg model to more exotic ones, like the Kitaev 
model 0. In designing effective spin systems different 
tools like polar molecules Q and tilted optical lattices Q 
can also be used. In recent years there has been a lot 
of interest in creating artificial Abelian and non-Abelian 
gauge fields in cold-atom systems [l(| and successful ex- 
perimental realizations of synthetic magnetic [ll L 
electric field [l3| and spin-orbit coupling (SOC) 
been reported. The role of SOC in cold atoms has been 
extensively studied following the theoretical proposal for 
the creation of artificial SOC [l|l [Tit and rich phase di- 
agrams have been found in BECs [17H19| and fermionic 
systems [23422 1 - 




In this letter we combine optical lattice and SOC 
which, in the deep lattice regime, leads to tight-binding 
description with non-zero "spin-flip" hopping between 
neighboring sites [23[. We show that in the Mott- 
insulator phase with integer filling the system is de- 
scribed by an interesting effective spin Hamiltonian 
which is a combination of Heisenberg model, quantum 



compass model and Dzyaloshinskii-Moriya terms. We 
note that combination of an optical lattice and SOC 
has already been considered with a purpose of study- 
ing superfluid-insulator transition [2JL topological phase 
transitions [25[ and BEC dynamics |26|. In the context 
of solid-state physics spin models resulting from Mott- 
insulator s with strong SOC were studied in Ref. [27H30| . 

We study a two-dimensional system of pseudospin-1/2 
atoms on a square optical lattice with synthetic SOC. 
The single-particle physics is described by the Hamilto- 



H a = 



£-+V x sin 2 {Kx) + V y sin 2 (Ky) 
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with p and m being the atomic momentum and mass; 
V x , V y the lattice depth in x and y direction, K = n/a 
(a is the lattice spacing), 1 the 2x2 unit matrix, &{ 
Pauli matrices; a and j3 characterize the SOC. Since 
([T]) is invariant under lattice translations, its eigenstates 
have Bloch-wave form: ijjk,n(^) = exp(ik • r)uk,n(r) and 
5k, n (r + R) = 4,n(r), R being a lattice vector. We are 
interested in the deep lattice regime in which pairs of en- 
ergy bands are well separated and the low-energy physics 
is captured by the lowest pair of bands which touch at 
k = (0,0), (0,-ftf), (K, 0) and (K,K) in the energy spec- 
trum. In this regime the system is well described by the 
tight-binding approximation in which the Hilbert space is 
spanned by states localized on individual lattice sites and 
the tunneling exists only between nearest-neighbor sites. 
There are two localized states per site (IW 7 ^), IW- 2 ,)), 
hence we have two effective particle species. This is the 
most general tight-binding description of (JlJ: 

Ht = £ J2 &i T ?' j)a ^ + H - c -] . ( 2 ) 
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where a£ i (a r .i) creates (annihilates) a particle in the 



state \W*), T- 



-(W*\H \WL 
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arc the tunnel- 



ing matrices and rj^ = cry. In finding the elements of 
T 7 corresponding to Ho, we use Peierls substitution and 
"localized Wannier states (LWS) method". We write 
SOC in a gauge- field form: p 2 / (2m) + aa x p x + /3a y p y = 
(p — A) 2 /(2m) + const, with A = (—maa x ,—m(3a y ) and 
may use Peierls substitution to find tunneling matrices 



y = x,y 



(3) 



where i 7 are tunneling coefficients in the 7-direction in 
the absence of SOC, 9 X = ira'/2, 9 y = 7r/3'/2; a' and 
P' are dimensionless SOC strengths: a' = 2ma/(HK), 
j3' = 2mP/{hK). However, Peierls substitution is only an 
approximation, valid for SOC weak with respect to the 
kinetic plus lattice part of Hq. This is the case when a' <C 
2^VjE^, p < 2^VjE^, where E R = h 2 K 2 /(2m) is 
the lattice recoil energy. While the SOC is quite weak 
in solid-state systems, in cold-atom systems it is typi- 
cally very strong and in that case Peierls substitution is 
not completely valid. For example, if we combine opti- 
cal lattice with spacing a = 410 nm and Rashba SOC 
generated by a scheme described in Ref. [3l|, we obtain 
a' ~ 1, while the usual experimental values of lattice 
depth arc \/V x , y /En ~ 2 — 5 (for smaller values of V XiV 
the tight-binding approximation is not valid). For the 
SOC scheme experimentally realized |l4[ a' « 1.4. Since 
the validity condition for Peierls substitution is not com- 
pletely satisfied, we calculate tunneling matrices using 
LWS method which is more involved and requires nu- 
merical approach but it does not contain any approxima- 
tion. In a system with N x N y sites and periodic bound- 
ary conditions, Wannier states for a single band are de- 



fined as [32J, l33|: |W, 



where |i/>k) are Bloch states and ip (k) is an arbitrary 
phase. In the absence of SOC it is possible to obtain 
maximally LWS by varying the phase <p(k) of each Bloch 
state |f/>k) [33j and in the deep lattice regime these max- 
imally LWS are well localized on individual sites. In 
the presence of SOC it is generally not possible to have 
Wannier states localized on individual sites if these are 
constructed from Bloch states of a single band. There- 
fore we consider generalized Wannier states introduced 



in Ref. [H]: |W? 



^N*N y 

where M mn (k) are 2x2 unitary matrices which mix Bloch 
states of the two bands. We find maximally LWS by 
minimizing the functional Q = J2n=i W 2 )" - ( r >«] with 
respect to matrix elements of M m „(k) ((...)« is an ex- 
pectation value associated to \Wq )). The minimization 
is done numerically and example of an algorithm is given 
in Ref. Numerical results show that the tunneling 

matrices still have the form given in ([3]) but now the 
parameters t x , t y , 9 X and 9 y are some more general func- 
tions of V x , V y , a, /3, K and m. It can be shown that 
the structure of tunneling matrices @ follows from the 
symmetries of Ho ■ 
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FIG. 1: (color online) Tunneling amplitudes Ty ' and T, 
obtained by Peierls substitution (full red and dash-dotted 
blue line) and localized Wannier states method (full dots 
and empty circles) for different strengths of Rashba SOC and 
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Peierls substitution has the advantage to give an ana- 
lytical form for tunneling matrices, however it does not 
give any information about the Wannier states, whereas 
the LWS method explicitly gives states \ W}), \W 2 ) which 
is important in interpreting the experimental data. In 
Fig. [T] we compare tunneling amplitudes in the Rashba- 
coupling case obtained by Peierls substitution and LWS 
method for V x = V y = 10 Er. They show excellent ac- 
cord for small a' and sizable differences for larger ones. 

Cold atoms in optical lattices are described by the tun- 
neling Hamiltonian ([2]) plus interactions 0]: 



V 



(4) 
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where number of particles in state \W£), Uij are 

interaction coefficients and :(...): denotes normal order- 
ing of creation and annihilation operators. 

We are interested in the Mott-insulator regime with 
tx/Utj -C 1, ty/Uij <C 1 and integer number v of atoms 
per site (y = 1 for fermions and any integer v for bosons). 
In this case interactions (0| are the dominant part in the 
full Hamiltonian H = H^ + V and we may treat the prob- 
lem pcrturbatively by taking V as starting Hamiltonian 
and Ht as perturbation. The ground state of V is a state 
with uniform distribution of atoms and the ground state 
degeneracy is very large since there are two states per 
site that atoms can occupy. The perturbation ifp cou- 
ples the ground-state manifold and excited states of V 
and the resulting low-energy effective Hamiltonian can 
be calculated in the second order of perturbation the- 
ory 0,0: 



(H e s)aP = - E 



(Ht)c:j(Ht, 
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where a and /? label states in the ground-state manifold, 
while 7 labels the excited states of V. 
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Wc first calculate the effective low-energy Hamiltonian 
for fermions. Since two fermions cannot occupy the same 
quantum state the only interesting regime is when v = 1. 
The only relevant interaction coefficient is Ui 2 = U and 
the excited states of V are those with two fermions of 
different species at the same site. Now it is convenient to 

+ 4.2«r,l)/2, 

OA 
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introduce isospin operators: = (a], ia T j 

S y = — i(aj x a r ,2 — 2 a r ,i)/2 and S" — (n r 
Using ([5} we obtain 



r' — r— r; 

+ E 
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where J£ = Jicos(20i), J l cm = 2J i sm z 0i, D 2 
J x sin(2<y £, = J a sin(20 v )j7, J 4 = 4tf/U and 6>, are 
introduced in The Hamiltonian is a combination of 
Heisenberg model, compass model and Dzyaloshinskii- 
Moriya-type terms; for 9 X = 9 y = (no SOC) and 



Jx = Jy we recover the Heisenberg model 

Next we find the effective low-energy Hamiltonian for 
bosons, and now the number of atoms per site v can be 
greater than one. The calculation for any v and general 
interaction coefficients is very cumbersome, however it 



simplifies for SU mn <^ U, where U„ 



U +SU mn or when 



v = 1 [5| . Wc express the Hamiltonian in terms of spin- 
v /2 operators defined in the same way as in the previous 
case. For 8U mn *C U, the Hamiltonian in the first order 
of 8U/U is H vB = -FF+Er [ A {SvY-hvS*] where A = 
(U n + U22 - 2E7ia) /2 and h v = -(v-l) (U n - U 22 ) /2. 
For v = 1 and general E/ mn wc obtain 



= E 



t—x.y r' —r—T), 

+ D^-(S r xS r 



T'i q q _i_ 7* j_ 7* C 12 
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where J\ = (u x + u 2 - 2)J£, J* m = 2Jisin 2 i , = 
(ui + it 2 )D l /2, fti = (ui - u 2 ) (J x cos 2 9 X + J y cos 2 # y ), 
Mi = U12/U11, u 2 = Ui 2 /U 22 and J x , J y are given in © 
with {/ replaced by U\ 2 . Since the atoms usually used in 
experiments have almost spin-independent interactions 
we assume Uij = U (u\ = u 2 = 1): this simplifies the 
Hamiltonian and yields Hib = —Hp. 

We intend to find the classical zero-temperature phase 
diagram of ([7]) with u\ = u 2 = 1. Some previous papers 
presented models combining Heisenberg, Dzyaloshinskii- 
Moriya and compass-model interactions [27H29| but they 
did not provide a complete phase diagram, neither at a 
classical level, usually considering only small SOC. In our 
approach we treat the spins as classical vectors and we 
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FIG. 2: (color online) Spin textures: (a) ferromagnet ((^ = 
= 0); (b) spiral wave (6 X — 0.5, 8 y = 0.2); (c) vortex phase 
(0x = 0y = 1); (d) antivortex phase (0* = 2.14, y = 0.96); 



'L 



(e) stripes 



1.6, 8 y = 0.7). Phases (a), (b) and (e) are 



coplanar and shown in a two-dimensional representation. 



aim to find the configurations {S^} which minimize the 
energy, with constant J x = Jy- We did our computations 
usually on 60 x 60-site lattices and finite-size effects are 
negligible. In Fig. [2] we show the phases and in Fig. [3] 
the corresponding phase diagram. We obtain two Ising- 
type phases [ferromagnet (Fig. [2^) and stripes (Fig. [2};)], 
coplanar spirals (Fig. ^jp) and three-dimensional ordered 
phases with vortices (Fig. or antivortices (Fig. [2J1). 
In describing our results, it is helpful to focus on a so 
called "basic region" given by the triangle 9 y < 6 X < ir/2: 
the solutions for other parameters can be obtained by 
simple mappings, e.g. ground-state configurations in 
X < 0y < 7r/2 region are obtained by simultaneous 7r/2- 
rotation of spins and sites of ground states in the "basic 
region". Upon activating SOC, the ferromagnet is imme- 
diately replaced by spiral waves, whose spatial periodicity 
reduces from several to three sites upon increasing X and 
y ] we found both commensurate and incommensurate 
waves. When the compass-model term becomes domi- 
nant over the Dzyaloshinskii-Moriya one, another copla- 
nar phase appears, the ferromagnetic stripe order, either 
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FIG. 3: (color online) Classical phase diagram of the Hamil- 
tonian Hib (see text for details): ferromagnet (black corner 
dots); spiral waves [dark orange (commensurate with four- 
sites periodicity), light orange (commensurate with three-sites 
periodicity), red (others)]; stripes (yellow); vortex phase (dark 
blue) and antivortex phase (light blue). 



directly or via the three-dimensional ordered phases. We 
always find non degenerate classical ground states, ex- 
cept along the dashed lines (Fig. [3j) which indicate points 
in the parameter space with a continuous degeneracy of 
classical ground states. However, we expect this degen- 
eracy to be removed by slight deviations of the realistic 
engineered SOC with respect to the Rashba-Dresselhaus 
form of the coupling |31[. The dashed lines also repre- 
sent the boundaries between phases with stripes of differ- 
ent orientation, i. e. between the phase shown in Fig. &p 
and the one obtained by rotating the sites and spins of 
the latter by tt/2 around the z-axis. The vortex phase 
(Fig. |2t) takes place along the diagonal 6 X = 9 y : vor- 
tices are left-handed in the region with smaller SOC and 
right-handed in the one with larger SOC. The antivortex 
phase (Fig. [3Ji) is found along the diagonal 9 y = n — 9 X 
and the configuration (d) is obtained from the phase (c) 
by a transformation which reflects sites (but not spins) 
with respect to the x-axis. For a better identification of 
the phase properties, we consider their behavior with re- 
spect to the breaking of the translational symmetry of 
0. While all the phases (except the ferromagnet) break 
this symmetry, they do it in a different way, i.e. the 
stripe phase in FigfSh is not invariant under one-lattice- 
sitc translation along x-direction, but it is invariant un- 
der two-lattice-sites translations in x-direction and un- 
der one-lattice-site translation in y-direction; the phases 
with vortices or antivortices are not invariant under one- 
and two-lattice-sites translations in x and y-direction, 
but they are invariant under three-lattice-sites transla- 
tions. Then we can understand the evolution from stripe 
to vortex phase as a transition in which two-latticc-sites 



translational symmetry becomes broken. The same rea- 
soning applies for the rest of the phase diagram. 

It is important to emphasize that the classical analysis 
yields no gapless modes in the whole parameter region 
except for the diagonal lines (9 y = 9 Xl 9 y = tt — 9 X ) 
in the stripe region. The absence of these gapless modes 
provides stronger guidelines for further analysis in a semi- 
classical or quantum approach. 

In summary, we studied the effects of spin-orbit cou- 
pling in cold atoms in an optical lattice in the Mott- 
insulating regime. We derived the tight-binding model 
using Peierls substitution and Localized Wannier State 
method and obtained the effective low-energy Hamilto- 
nian for fermions and bosons: this takes the form of an 
exotic spin model with Heisenberg, compass-model and 
Dzyaloshinskii-Moriya interactions. We determined the 
classical phase diagram for this model and showed that 
the interplay between the different interactions is respon- 
sible for a large variety of phases: ferromagnet, spirals, 
stripes, three-dimensional vortex and antivortex phases. 
We expect that our classification of ground states could 
generally survive in a quantum approach; in fact, ex- 
cept for some particular cases we mentioned in the dis- 
cussion, on the classical level there are no degeneracies 
which would be lifted by quantum fluctuations. 
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